Functional analysis of the airways after pulmonary lobectomy through computational fluid dynamics

Pulmonary lobectomy, which consists of the partial or complete resection of a lung lobe, is the gold standard intervention for lung cancer removal. The removal of functional tissue during the surgery and the re-adaptation of the remaining thoracic structures decrease the patient's post-operative pulmonary function. Residual functionality is evaluated through pulmonary function tests, which account for the number of resected segments without considering local structural alterations and provide an average at-the-mouth estimation. Computational Fluid Dynamics (CFD) has been demonstrated to provide patient-specific, quantitative, and local information about airways airflow dynamics. A CFD investigation was performed on image-based airway trees reconstructed before and after the surgery for twelve patients who underwent lobectomy at different lobes. The geometrical alterations and the variations in fluid dynamics parameters and in lobar ventilation between the pre and post-operative conditions were evaluated. The post-operative function was estimated and compared with current clinical algorithms and with actual clinical data. The post-operative configuration revealed a high intersubject variability: regardless of the lobectomy site, an increment of global velocity, wall pressure, and wall shear stress was observed. Local flow disturbances also emerged at, and downstream of, the resection site. The analysis of lobar ventilation showed severe variations in the volume flow rate distribution, highlighting the compensatory effects in the contralateral lung with an increment of inflow. The estimation of post-operative function through CFD was comparable with the current clinical algorithm and the actual spirometric measurements. The results confirmed that CFD could provide additional information to support the current clinical approaches both in the operability assessment and in the prescription of personalized respiratory rehabilitation.

Mesh selection. The convergence analysis showed no relevant differences between the three meshes as reported in Fig. 1c for a representative subject in the pre and post-operative models. Consequently, the meshes with the lowest number of elements were adopted. The average number of elements of the mesh was 3.2 × 10 6 (± 0.5 × 10 6 ) and 3.6 × 10 6 (± 1.1 × 10 6 ) for the pre and post-operative models, respectively. The average skewness was equal to 0.87 ± 0.02. Table 1. Study population. For each patient, the location of the lobectomy site, sex, age and clinically measured forced expiratory volume in the 1st second (FEV1) before (preFEV1) and after (postFEV1) the surgery is reported. the trachea and the main bronchi and the cross-sectional area (CSA) of the remaining lobar segment adjacent to the resection site measured in the middle point of the segment and perpendicular to the centerline are reported in Table 2 for the pre and post-operative models, divided according to the different sites of lobectomy. We observed specific trends in upper lobectomies (both RUL and LUL): the angle between the trachea and the main bronchi on the operated side tends to decrease while showing a mild increase on the contralateral side. Increments of θ R were also observed in RLL lobectomies with an average increase of 10.94 (7.72-13.54) %, while θ L showed mild variations. None of the differences in θ R and θ L between the pre and post-operative models were statistically significant. For both the RUL and LUL lobectomy, the secondary bronchi are twisted upward and severely bent downstream of the suture zone, as can observed in Fig. 2a and Fig. 2c. The CSA of the adjacent segment in the post-operative models resulted significantly lower than the pre-operative ones (p = 0.01) for all the subtypes of lobectomies considered. In general, higher reductions were observed for upper lobectomies: the median decrement was equal to 21 (18)(19)(20)(21)(22)(23)(24)(25)(26)(27)(28) % and 34 (26)(27)(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38)(39)(40)(41)(42)(43) % for RUL and LUL lobectomies, respectively. Excluding patient 11, lower lobe lobectomies showed a cross-sectional reduction of 0.50 (0.40-3.60) %. For subject 11, marked cross-sectional variations were observed with a value equal to -86%. Results concerning this subject are detailed in the Supplementary Material Section. In LLL and RLL lobectomy cases, the geometrical alterations of the remaining lobes appeared to be less marked. No evidence of bronchial stenosis was observed in the remaining lobar segments. Moderate deformations of the RML and of the lingula were observed for RLL and LLL lobectomies cases, respectively.

Fluid dynamics parameters.
Globally and regardless of the resected lobe, the post-operative configurations showed, given the same inlet flow, an increment in the magnitude of velocity, wall pressure, and wall shear stress for all 12 patients. The peak values of the variables of interest are reported in Table 3. In the comparison between the fluid dynamics parameters in the pre and post-operative models, statistically, significant differences were observed between the peak values of pressure (p < 0.001), velocity (p = 0.0017), and wall shear stress    Global peaks of wall shear stresses generally occurred at the first two generations of the tree, although for LUL and RUL lobectomies, comparable values were found in the first bronchial generation downstream of the suture. Finally, also the pressure drop between the inlet and the outlets significantly increased (p < 0.001) in the post-operative models. Significant negative correlations were identified between the postoperative CSA of the remaining segment and the maximum values of pressure (r = -0.881, p < 0.001), velocity (r = 0.860, p = 0.004), wall shear stress (r = 0.697, p = 0.012). The lower the CSA, the higher the pressure drop (r = -0.888, p < 0.001). No correlation was identified between the values of the θ R and θ L and the fluid dynamic parameters. In Fig. 3, the local comparison of the fluid dynamics on the surgery site for four patients representative of each site of lobectomy is considered. Local maxima of pressure, velocity, and wall shear stress were present at the resection site level and in correspondence to the remaining lobar segments for all the patients under analysis. These local maxima were more marked for those patients who underwent an upper lobe lobectomy and presented evidence of cross-sectional reduction at the levels of the remaining secondary bronchi and severe rearrangements of the remaining lung structures as previously described. For RLL and LLL lobectomy cases, local increment of pressure, wall shear, and velocity were still observable despite the absence of marked geometrical alteration due to the redistribution of flow in the upper lobes. Interestingly, consequences of the lobectomy were also observable on the contralateral side of the tree, where no surgical intervention was performed. Higher velocity, pressure, and wall shear stress values were found independently from the lobectomy site. In one patient (#5, reported in the supplementary material), whose trachea assumed a deformed configuration after surgery, local maxima were also present. No statistically significant differences in terms of percentual variation of pressure, velocity, or wall shear stress were observed for the patients under analysis if grouped for the affected lung (right/left) or position (upper/lower). Figure 4 shows the distribution of VFR for the subjects under analysis grouped according to the site of lobectomy. In the pre-operative models (white bars), the right lung received a higher amount of flow. Lower lobes showed higher flow than the upper ones, and the RML was the lobe constantly receiving less flow. The post-operative models revealed considerable rearrangements of the VFR distribution for all the typologies of lobectomy considered. In the comparison of the overall population, significant differences arose between the pre and post-operative models for both the operated lung (p < 0.001) and the contralateral lung (p < 0.001). However, if the lobe to be removed is neglected in calculating the VFR to the operated lungs, the difference between the pre and post-operative configurations is not significant (p = 0.121). Finally, in the post- www.nature.com/scientificreports/ operative models, a statistically significant imbalance was found between the operated and the contralateral lung with a significant increment (p < 0.001), in the contralateral lung. This was not observed in the pre-operative comparison between the lung affected by the tumor and the contralateral one (p = 0.384). Table 4 summarizes the value of the clinically measured FEV1 and the estimated post-operative function with both the anatomical formula and the CFD approach. The FEV1 estimation through CFD (Eq. (2.4)) showed a significant correlation with the reference clinical estimation through Eq. (2.3) (r = 0.95, p < 0.001), and with the actual clinically measured post-operative FEV1 (r = 0.806, p = 0.002). For the subjects under analysis, both the CFD method and the anatomical formula tend to over-predict the clinically measured FEV1 after the surgery.

Discussion
In the present work, CFD was applied to evaluate the post-operative fluid dynamic alterations following pulmonary lobectomies performed in different lobes. Patients' specific geometries were reconstructed from CT data, the airflow simulated, and the results analyzed to identify fluid dynamics parameters and VFR variations. Computational results were used to estimate the post-operative function (postFEV1CFD), which was compared against the current clinical standard algorithm (postFEV1anatomic) and the actual post-operative measurements (postFEV1). The present study allowed for a more quantitative assessment of the pre-operative condition and an enhanced understanding of the post-operative alterations. The information of pre-operative flow rate can be used, according to Eq. (2.4), to estimate the post-operative function and provide additional information concerning patient operability without the actual need for post-operative data. In contrast to the post-operative estimation performed with the anatomical formula, CFD accounts for local fluid dynamics characteristic of each branch that will eventually be resected during the surgery.
Regarding the numerical simulations' settings, it was chosen to simulate a steady-state peak inspiratory condition. The choice of a steady-state simulation was a reasonable simplification considering the results related to   25,26 . A uniform pressure was set at the outlet boundaries not to constrain the airflow partition and the lobar flow rates 15 . The k-ω SST model has been widely adopted for airflow simulation in the airways, and it has been demonstrated to provide an optimal trade-off between the accuracy of the prediction of low turbulence flow in the respiratory system and the availability of computational resources 15,23,27 . The choice of considering the air as a Newtonian fluid at ambient conditions is motivated by the fact that the Mach number (the ratio between the velocity of air and the speed of sound) in the airways is approximately 0.3 and that temperature and humidity variations within the physiological range of the respiratory system cause negligible variations on the air properties with respect to environmental    [29][30][31][32] , can be considered acceptable due to the relatively high percentage of cartilage in the airway tissue of the central airways 33 .
In terms of geometrical characterization of the pre and post-operative models, we observed severe alterations of the remaining pulmonary structures in the case of upper lobectomies. Although present, the structural modifications were less relevant for lower lobe lobectomies and mostly limited to the RML and the lingula shift to the vacant area, in accordance with previous studies 34 . The alterations were also relevant at the trachea and main bronchi level in two subjects, as described by previous studies 24,35 . The significant decrement of the cross-section downstream the surgery site is confirmed by studies for the LUL cases 23,24 , and the results of the proposed work suggest that this consideration is also verified in the RUL lobectomies cases. Lower lobe lobectomies, on average, showed lower cross-section variations with the only exception of subject 11. As reported in the supplementary material section, the matching with the clinical course revealed that the patient had suffered from pneumonia with pleural effusion. Our analysis readily detected this abnormal clinical course, and the results were widely different from other cases. Indeed, if this subject is neglected, the median reduction of CSA in the case of lower lobectomies is equal to 0.60 (0. 40-5.53). In contrast to Gu et al. 23 , no systematic increment of θ R and decrement of θ L was observed in the LUL post-operative cases, possibly because of the limited size of our population. Although the characterization and the role of structural modifications on post-operative pulmonary function is still an open field of research 36 , the pre and post-operative values measured angles in our study are in line with the expected intersubject variability, also accounting for the sex-related differences between women and men 37 .
The numerical simulation results showed that, under the same inlet flow rate, increased mechanical stimuli emerged in the post-operative models with increased velocity, pressure, and wall shear stress in the whole bronchial tree. These fluid dynamic alterations are derived from the resection of an entire lung lobe and the geometrical adaptation of the remaining structures, especially for upper lobe lobectomies. Limited to LUL lobectomies, our results are consistent with the two previous studies characterizing this type of lobectomy 23,24 . With respect to Gu et al. 23 , relatively higher variations of fluid dynamic parameters are observed in our cohort of patients that underwent LUL lobectomy, although this can be possibly ascribed to the higher inlet flow rate considered in the simulations of the presented work. The obtained values are indeed in line with Tullio et al. 24 , who considered the same inlet flow of the presented study. For upper lobes lobectomies, the increment of pressure, velocity, and wall shear stress can be associated with the severe bronchial stenosis observed and the general re-adaptation of the remaining lower lobes. For lower lobes lobectomies, despite the presence of mild stenosis of the bronchi upward the suture site, it seems more likely that the fluid dynamics variations are associated with the redistribution of the flow that physiologically is much higher the lower lobes than to the upper ones. The increment in pressure drop observed after the surgery, independently from the lobectomy site, can suggest a higher fluid dynamic resistance to airflow and a consequent pressure-driven increase during inspiration 18,38,39 . The result is consistent with previous computational studies on lobectomy 23,24 . The quantitative and local information on fluid dynamics parameters provides a better understanding of the alteration induced by the re-adaptation of remaining lung structures and has a physiological and clinical impact. Indeed, the increment in the values of wall shear stress caused by the increased local velocity at the level of the surgical site can impact the endothelial cells of the bronchial tissue 40 , and represent a risk especially in case of higher inlet flow rates, for example during exercise conditions.
In terms of flow distribution, the results in the pre-operative models are consistent with previous airflow dynamics and particle deposition studies. The right lung receives more flow due to its bigger size if compared to the left one 41,42 , although, as expected, this difference is not statistically relevant. The lower lobes receive a greater portion of inhaled flow rate with respect to the upper lobes, while the RML, being the smallest lung lobe, receives the lowest amount of flow rate 15,17,43 . In the post-operative models, a significant increment in the contralateral lung flow with respect to the operated lobe was observed, as confirmed by previous computational studies limited to the LUL lobectomies cases 23,24 . The presented findings suggest that this trend is also confirmed for RUL, LLL, and RLL resections. However, contrary to Gu et al. 23 , the remaining structure in the operated lung does not seem to undergo statistically significant variations in the inflow, as in LUL lobectomies. This supports the idea that the contralateral lung tends to undergo the major flow variations and takes most of the flow physiologically meant for the resected lobe. These findings are supported by clinical studies on pulmonary function after lobectomy, which evidenced a reduction in volume variations and regional ventilation in the remaining lobes of the operated lung 44,45 .
Morphological re-adaptations and the mediastinum shifting, make the measurement of lobar ventilation distribution through functional imaging particularly challenging in pulmonary lobectomy surgery 44 . CFD can therefore provide clinically relevant support to address this issue. In addition, the possibility to obtain patientspecific information on fluid dynamics and lobar ventilation could become a valuable asset in the prescription of personalized respiratory physiotherapy for better post-operative management.
The comparison presented in Table 4 among the estimations of the post operative function with the CFD (postFEV1CFD), the anatomical method (postFEV1anatomic) and the clinical measured values (postFEV) showed an overall good consistency. However, it appears that the CFD approach does not currently provide a better estimation of the actual clinical measure. This can be possibly related to the relatively limited number of subjects under analysis and to relatively simplified settings adopted in the numerical simulations, in particular in terms of inlet and outlet boundary conditions.
The proposed work is affected by some limitations deriving from the available input data and the computational resources at disposal. Several studies have demonstrated the efficacy of adopting patient-specific lobar flow rate of the patient as outlet boundary condition 32,[46][47][48] : this requires functional imaging or at least the evaluation of volume variation between inspiratory or expiratory scans. Due to the retrospective nature of the collected data, no expiratory scans were available, as the current clinical follow-ups for these patients include only inspiratory CT images. For this reason, as most authors in literature 15 www.nature.com/scientificreports/ in the simulations. As an alternative, impedance-based boundary conditions that account for the impedance of the parts of the airway tree that cannot be resolved (either for computational reasons or for intrinsic limitations of the CT), can be applied to provide physiological boundary conditions for flow in the tracheobronchial region 52 . Besides, to better characterize the fluid dynamics in the proximity of the trachea inlet and prevent artifacts related to the application of the boundary condition, it would be beneficial to include flow extension to allow for flow development or directly imposing developed profile at the inlet. No direct validation of the results related to CFD was performed on these patients. Nevertheless, the proposed numerical approaches have been validated and have been demonstrated to reasonably agree with functional measurements from imaging (e.g., in Luo et al. 15 ) and to more complex models such as DNS 26 . In addition, although consistent with previous literature, the assumption of a steady-state solver neglects possible unsteadiness. Finally, in our analyses, only the limit case peak inspiratory flow under normal breathing was considered because it was expected to provide the highest mechanical stimuli in the airway walls 25,26 . Future directions of development will include the adoption of a dynamic inlet flow rate mimicking a physiological-like breathing waveform 48 and the investigation of different breathing conditions during exercise or mechanical ventilation, including the actual breathing waveform of the patients or flow rate data from literature 45 . Unsteady simulations can also be performed at the same inlet flow rate to account for unsteadiness. Prospective studies on a selected cohort of patients, including expiratory scans and functional imaging can be designed both to apply patientspecific boundary conditions and to validate the numerical solution on the specific dataset. In addition, it would be possible to perform the CFD analysis of the expiratory phase. Ideally, the acquisition protocol should also include the patient's upper airways for a more realistic representation. Also, to better understand the post-surgery remodeling, the geometrical characterization of the airways in the pre and post-operative conditions can be further developed with more comprehensive metrics such as the regional volumetric changes and the evaluation of different cross-sections along the bronchial tree. Finally, an extended cohort of subjects, including patients who underwent middle lobe lobectomies and bilobetomies, should be analyzed. Finally, an extended cohort of subjects, including patients who underwent middle lobe lobectomies and bilobetomies, should be analyzed. The inclusion of more patients, together with the evaluation of more complex numerical modelling approaches, would also be highly beneficial in the assessment of the degree of accuracy of post-operative function estimation through the CFD approach compared to the current clinical algorithm.

Conclusions
A CFD approach to evaluate the consequences of post-operative anatomical alterations on airflow dynamics in the airways model of patients that underwent pulmonary lobectomy of different lobes has been proposed. The airflow disturbances were quantified both globally and locally. To the best of our knowledge, this is the first study that applies CFD to subjects who underwent lobectomies in different sites. The results highlighted that alterations related not only to the resection of an entire lobar structure but also to the subsequent re-adaptation of the remaining lung structures. The regional flow rate to the lobes was computed, suggesting that compensatory effects occur at the level of the contralateral lung rather than in the remaining operated lung lobe(s). The estimation of the post operative function from the simulations (postFEV1CFD) was comparable with the results of the current clinical algorithm (anatomical calculation, postFEV1anatomic) and the actual clinical data (postFEV1). This approach will provide additional information concerning patient operability, accounting for actual patientspecific morphology and functionality alterations and for local fluid dynamics characteristic of each branch that will eventually be resected during the surgery.

Methods
Study design. This retrospective study was approved by the Fondazione IRCCS Ca' Granda -Ospedale Maggiore Policlinico of Milan institutional review board (approval number: 2583), and informed consent was obtained from all patients. All methods were carried out in accordance with relevant guidelines and regulations. A query of the institution's radiology database for chest CT images acquired for clinical purpose in breath-hold at full inspiration one month before and six months after surgery was performed for early-stage lung cancer patients, undergoing pulmonary lobectomy, aged > 18 years old. A subgroup of CT examinations, which presented technical acquisition parameters suitable for the analysis, were selected for the study. All identifying information was removed from the CT images before the analysis. Clinical details related to the transplantation and spirometry results were obtained through the electronic medical record. Exclusion criteria were middle lobectomy, bilobectomy, CT performed more than one month before or more than six months after surgery, patients with no clinical details available both before and after surgery. Chest CT images acquired at Thoracic Surgery and Lung Transplantation Unit of Fondazione IRCCS Ca' Granda -Ospedale Maggiore Policlinico of Milan, Italy.All methods were carried out in accordance with relevant guidelines and regulations. The following CT scan (SOMATOM Definition dual-source CT; Siemens, Forchheim, Germany) settings were adopted during image acquisition: tube voltage = 100-120 kV; tube current = 142-617 mA according to body mass index, matrix size = 512 × 512, slice thickness = 0.62-1 mm, in-plane resolution = 0.62-0.90 mm; axial resolution = 0.60-1.0 mm. Reconstruction was obtained with soft kernels (B30f., B31f.).
Spirometry, including forced vital capacity (FVC) and forced expiratory volume in 1 s (FEV1), were collected before and after the surgery, at the time of CT, and reported in terms of percentage predicted of normal values. The surgery was carried out in the period from January 2015 to June 2019. We divided the study population into 4 groups, one for each type of lobectomy: right upper lobectomy (RUL), right lower lobectomy (RLL), left upper lobectomy (LUL), left lower lobectomy (LLL). We consecutively selected 12 subjects meeting the inclusion/ exclusion criteria, up to reach 3 complete cases in each of the groups. www.nature.com/scientificreports/ Generation of the airway model. The 3D subject-specific airway tree models for the pre-operative and post-operative configurations were semi-automatically segmented in Mimics (Materialise NV, Belgium). The geometrical reconstructions were processed with 3-Matic (Materialise NV, Belgium) and ANSYS SpaceClaim (ANSYS Inc., Pennsylvania, USA) and export to ANSYS Fluent for meshing. The average volume of the reconstructed geometries was equal to 5.97 × 10 4 (± 2 × 10 4 ) mm 3 and 5.43 × 10 4 (± 1.99 × 10 4 ) mm 3 for the pre and post operative geometries respectively.

Mesh generation and convergence analysis.
A tetrahedral unstructured grid type was chosen due to the geometrical complexity. To resolve fluid properties in the region proximal to the wall, 15 inflation layers were created. The first layer thickness was set equal to 0.01 mm, corresponding to an estimated y + < 1. The quality of the generated mesh was evaluated through the maximum skewness criterion (equilateral volume skewness < 0.90). Three mesh refinements (M1, M2, M3) were generated for each geometry. The total number of elements for each mesh varied according to the complexity of the patient geometries. Over the whole population, for the pre-operative cases, the average number of elements was equal to 3.2 × 10 6 (± 0.5 × 10 6 ), 5 × 10 6 (± 1.1 × 10 6 ), 9.3 × 10 6 (± 2.4 × 10 6 ) for M1, M2, M3, respectively. For the post-operative cases, the average number of elements was equal to 3.6 × 10 6 (± 1.1 × 10 6 ), 5.7 × 10 6 (± 2.1 × 10 6 ), 10.0 × 10 6 (± 3.9 × 10 6 ) for M1, M2, M3, respectively. The convergence analysis was conducted considering the relative variations between the meshes regarding volume flow rate (VFR) at the outlets, with a tolerance inferior to 0.5%. Moreover, the local velocity profile was evaluated for the three meshes at a cross-sectional plane at 2 cm distance from the carina (Fig. 1a). Two axes, mutually rotated at 90° from each other, were chosen to plot the profiles as shown in Fig. 1b. Airflow simulation. The influence of oscillatory behavior on the flow was preliminarily investigated. The effects of unsteadiness on the average flow characteristics, computed through a steady solution, are modest if the Womersley (Wo) and the Strouhal (S) numbers are lower than 10 and 1, respectively 41 . The Womersley number, i.e., the ratio of unsteady forces to viscous forces, is calculated as: The Strouhal number, i.e., the ratio of unsteady forces to inertial forces, is calculated as: where f is the breathing frequency, υ is the air kinematic viscosity, D is the hydraulic diameter of the trachea, and u is the average velocity at the corresponding diameter. In the cases under analysis, f = 15 breaths per minute and υ = 1.79 × 10 -5 kg/ms were assumed. D was measured at the inlet cross-section for all the pre and post-operative models, while the velocity at the trachea was obtained by dividing the inlet flow rate by the inlet cross-section. The inspiratory flow was assumed to equal 0.5 L/s 53 , simulating the peak inspiratory flow during quiet breathing. Under those conditions, the average Womersley number was equal to 2.45 (± 0.35), and the average Strouhal number was equal to 0.012 (± 0.005). The Womersley and Strouhal numbers were equal to 2.38 (± 0.37) and 0.011 (± 0.005) for the post-operative cases. Consequently, steady-state simulations were performed. The finite volumes method was adopted for the numerical solution. Since the average Reynolds number was equal to 2606 ± 394 and 2707 ± 457 for the pre and post-operative cases, respectively, turbulence was modeled using the Shear Stress Transport (SST) k-ω model with Low Reynolds Number (LRN) correction with a turbulent intensity of 5% and a viscosity ratio (µ T /µ) of 10. A steady inspiratory flow corresponding to 0.5 L/s was set as a boundary condition at the trachea inlet. A uniform reference pressure was set at each outlet boundary, and a no-slip condition was imposed at the solid walls, considered rigid, stationary, and smooth. The air was modeled as an incompressible Newtonian fluid, with density and dynamic viscosity equal to 1.225 kg/m 3 and 1.79 × 10 −5 kg/ms. This is an acceptable assumption for the airflow in the human airways, where the Mach number (the ratio between the velocity of air and the speed of sound) is close to 0.3, and the variation air properties are negligible 28,41 . The hypothesis of rigid walls with the no-slip condition is well documented in the literature [29][30][31][32]41 and can be considered reasonably valid due to the relatively high percentage of cartilage in the airway tissue of central airways 33 .
The following settings were adopted: pressure-based solver; a second-order upwind scheme for the momentum equations; Green-Gauss cell-based method for gradients evaluation; and the semi-implicit method for pressure linked equations (SIMPLE) algorithm for pressure-velocity coupling. Convergence was set as residuals less than 10 −6 . Calculations stopped when the residuals converged, and the solution was evaluated to be stable. The results were reported in terms of velocity magnitude, wall pressure, and wall shear stress. The pressure drop was calculated as the difference between the mean pressure at the inlet of the trachea and the mean pressure of the outlets. Finally, the regional ventilation of each lobe was computed for both the pre-operative and post-operative cases by summing the VFR of the corresponding outlets.
Geometrical characterization. The right (θ R ) and left (θ L ) in-plane bifurcations angle between the trachea and the main bronchi were manually measured using the bifurcating nodes of the centreline reference points. The operation was performed on both the pre and post-operative geometries. Additionally, the cross-sectional area (CSA) of the remaining lobar segment adjacent to the resection site was measured in the middle point of the segment. The cross-section perpendicular to the centerline of the airway in the given point was selected. (postFEV1 anatomic ) has been performed in accordance with the current clinical guidelines 9,54 , using the anatomical, segment counting method 10,11 : where y is the number of functional or unobstructed lung segments to be removed, and z is the total number of functional segments. The total number of segments for both lungs is assumed to be equal to 10 (3 RUL, 3 in the RML, and 4 in the RLL) for the right lung and 9 to the left lung (5 in the LUL, 4 in the LLL). The estimation of the post operative funtion through the CFD (postFEV1CFD) approach was computed by multiplying the pre-operative function by the remained functionality of the lung estimated as the percentual flow to the remaining lobes, that will not be resected during the surgery (VFR remaning lobes ).
The post-operative function estimated with the anatomical formula (postFEV1anatomic) and the CFD method (postFEV1CFD) was compared against the actual clinical FEV1 (postFEV1) measured from the patients after the surgery.
Statistical analysis. Study population and mesh element number data are presented as mean ± standard deviation. Values of velocity, wall pressure, walls shear stress, and pressure drop are reported as median (25th percentile-75th percentile) unless stated otherwise. Relative variations were computed with respect to the preoperative condition as reported in Eq. (2.5): where X is the parameter of interest.
The assumption of normality was tested through Shapiro-Wilk Test. If normality was verified, paired t-test was adopted to confront data. For not parametric data, Wilcoxon matched pairs test was adopted. Pearson Correlation was used to evaluate correlations. The level of statistical significance was set at P < 0.05.